{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "59c0fbea-9eda-4477-be6c-32647cef1ec7",
   "metadata": {},
   "source": [
    "# CI Simulation (Supplementary Material)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12618aaa-7980-4ed3-8cd8-cd7ce3d3809a",
   "metadata": {},
   "source": [
    "- This notebook supplements the [confidence-intervals-for-ml.ipynb](confidence-intervals-for-ml.ipynb) with a case study.\n",
    "\n",
    "- Here, we are interested in seeing whether the true model accuracy (generalization accuracy) is actually contained in the confidence intervals.\n",
    "\n",
    "- For this, we create a synthetic dataset consiting of 10 million and 2 thousand data points for classification as shown in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8412450-153d-4788-89ed-0e5caa6e7b7b",
   "metadata": {},
   "source": [
    "## Large Synthetic Training and Test Sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91d86525-404f-4c82-80b5-38eafc26955f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_classification\n",
    "\n",
    "X, y = make_classification(\n",
    "    n_samples=10_002_000,\n",
    "    n_features=5,\n",
    "    n_redundant=2,\n",
    "    n_classes=2,\n",
    "    n_clusters_per_class=1,\n",
    "    random_state=123,\n",
    "    flip_y=0.25,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fd6fb5cc-24d1-47bd-8682-f4fd25f1eda4",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = X[:1_000]\n",
    "y_train = y[:1_000]\n",
    "\n",
    "X_test = X[1_000:2_000]\n",
    "y_test = y[1_000:2_000]\n",
    "\n",
    "X_huge_test = X[2_000:]\n",
    "y_huge_test = y[2_000:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c018df3-637b-41bd-b270-4935527fff94",
   "metadata": {},
   "source": [
    "- Note that the 1000 data points are used for training, the second 1000 data points are used for testing, and the remaining 10,000,000 data points represent the dataset we use to calculate the true performance of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a391f1e1-d3ab-45fe-81ee-bf5a48600afd",
   "metadata": {},
   "source": [
    "## True Generalization Performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7ce712bf-e9ba-4764-a3af-e03e1261939b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8472259"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "clf = DecisionTreeClassifier(random_state=123, max_depth=3)\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "acc_test_true = clf.score(X_huge_test, y_huge_test)\n",
    "acc_test_true"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68a4eee6-2422-4b9d-9d89-a7ffef36d730",
   "metadata": {},
   "source": [
    "## 1) Normal Approximation Interval Based on the Test Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "98aad95f-e372-4766-9912-5828f2f83720",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8289298133898874 0.8730701866101126\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats\n",
    "\n",
    "confidence = 0.95  # Change to your desired confidence level\n",
    "z_value = scipy.stats.norm.ppf((1 + confidence) / 2.0)\n",
    "\n",
    "\n",
    "clf.fit(X_train, y_train)\n",
    "\n",
    "acc_test = clf.score(X_test, y_test)\n",
    "ci_length = z_value * np.sqrt((acc_test * (1 - acc_test)) / y_test.shape[0])\n",
    "\n",
    "ci_lower = acc_test - ci_length\n",
    "ci_upper = acc_test + ci_length\n",
    "\n",
    "print(ci_lower, ci_upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9a7f37e8-2bf5-447f-be9c-9306c6881a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {\n",
    "    \"Method 1: Normal approximation\": {\n",
    "        \"Test accuracy\": acc_test,\n",
    "        \"Lower 95% CI\": ci_lower,\n",
    "        \"Upper 95% CI\": ci_upper,\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3e93a0b-db9a-4b84-ad5b-3e230d4e3a38",
   "metadata": {},
   "source": [
    "### 2) Out-of-Bag (OOB) Bootstrap; Bootstrapping Training Sets -- Setup Step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3cdd76d5-d7be-4fe4-9f08-29e81e6ae3d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "\n",
    "\n",
    "for i in range(bootstrap_rounds):\n",
    "\n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "\n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "\n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    acc = clf.score(boot_valid_X, boot_valid_y)\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "\n",
    "\n",
    "confidence = 0.95  # Change to your desired confidence level\n",
    "t_value = scipy.stats.t.ppf((1 + confidence) / 2.0, df=bootstrap_rounds - 1)\n",
    "\n",
    "\n",
    "se = 0.0\n",
    "for acc in bootstrap_train_accuracies:\n",
    "    se += (acc - bootstrap_train_mean) ** 2\n",
    "se = np.sqrt((1.0 / (bootstrap_rounds - 1)) * se)\n",
    "\n",
    "ci_length = t_value * se\n",
    "\n",
    "ci_lower_21 = bootstrap_train_mean - ci_length\n",
    "ci_upper_21 = bootstrap_train_mean + ci_length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4d8e2ba4-c453-4978-a215-55a6fcabfc6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.1: Bootstrap, 1-sample CI\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower_21,\n",
    "    \"Upper 95% CI\": ci_upper_21,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0959f31-e2a0-4aca-818b-0e7315ec2faf",
   "metadata": {},
   "source": [
    "### 2.2) Bootstrap Percentile Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "88c69d47-1cca-42e0-ac10-4e069b678be0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_lower_22 = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper_22 = np.percentile(bootstrap_train_accuracies, 97.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "18c7fcbd-958b-4942-9cef-fe6aa7dbcdb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "results[\"Method 2.2: Bootstrap, percentile\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower_22,\n",
    "    \"Upper 95% CI\": ci_upper_22,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76172a2c-c24e-4fba-8653-fa6da4d45b43",
   "metadata": {},
   "source": [
    "### 2.3) .632 Bootstrap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "946f60f0-b2af-43fb-9a56-34b46786a494",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "weight = 0.632\n",
    "\n",
    "for i in range(bootstrap_rounds):\n",
    "\n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "\n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "\n",
    "    clf.fit(boot_train_X, boot_train_y)\n",
    "    train_acc = clf.score(X_train, y_train)\n",
    "    valid_acc = clf.score(boot_valid_X, boot_valid_y)\n",
    "    acc = weight * train_acc + (1.0 - weight) * valid_acc\n",
    "\n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "da458747-9f9e-47ee-a410-261e0c69ef01",
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_lower_23 = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper_23 = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "\n",
    "results[\"Method 2.3: Bootstrap, .632\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower_23,\n",
    "    \"Upper 95% CI\": ci_upper_23,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d4d29b-6e4b-42db-b50f-6b8d6750c87a",
   "metadata": {},
   "source": [
    "### 2.4) .632+ Bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3483410-c40c-4efe-aa63-d7a7467659c4",
   "metadata": {},
   "source": [
    "- Unfortunately, this method is too computationally expensive for this dataset on a regular computer, which is why we skip it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0ac8e342-0d99-4c1f-ad4b-5217748f158e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nfrom itertools import product\\nfrom sklearn.base import clone\\nfrom sklearn.metrics import accuracy_score\\n\\n\\ndef no_information_rate(targets, predictions, loss_fn):\\n    combinations = np.array(list(product(targets, predictions)))\\n    return loss_fn(combinations[:, 0], combinations[:, 1])\\n\\nrng = np.random.RandomState(seed=12345)\\nidx = np.arange(y_train.shape[0])\\n\\nbootstrap_train_accuracies = []\\nbootstrap_rounds = 200\\nweight = 0.632\\n\\n\\ncloned_clf = clone(clf)\\nfor i in range(bootstrap_rounds):\\n    \\n    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\\n    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\\n    \\n    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\\n    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\\n    \\n    cloned_clf.fit(boot_train_X, boot_train_y) \\n    train_acc = clf.score(X_train, y_train)\\n    valid_acc = cloned_clf.score(boot_valid_X, boot_valid_y)\\n    \\n    gamma = (no_information_rate(\\n        y, cloned_clf.predict(X),\\n        accuracy_score))\\n    R = (valid_acc - train_acc) / (\\n        gamma - train_acc)\\n\\n    weight = 0.632 / (1 - 0.368*R)\\n    \\n    \\n    acc = (weight*train_acc + (1. - weight)*valid_acc)\\n    \\n    bootstrap_train_accuracies.append(acc)\\n\\nbootstrap_train_mean = np.mean(bootstrap_train_accuracies)\\n\\n'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "from itertools import product\n",
    "from sklearn.base import clone\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "\n",
    "def no_information_rate(targets, predictions, loss_fn):\n",
    "    combinations = np.array(list(product(targets, predictions)))\n",
    "    return loss_fn(combinations[:, 0], combinations[:, 1])\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_train.shape[0])\n",
    "\n",
    "bootstrap_train_accuracies = []\n",
    "bootstrap_rounds = 200\n",
    "weight = 0.632\n",
    "\n",
    "\n",
    "cloned_clf = clone(clf)\n",
    "for i in range(bootstrap_rounds):\n",
    "    \n",
    "    train_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    valid_idx = np.setdiff1d(idx, train_idx, assume_unique=False)\n",
    "    \n",
    "    boot_train_X, boot_train_y = X_train[train_idx], y_train[train_idx]\n",
    "    boot_valid_X, boot_valid_y = X_train[valid_idx], y_train[valid_idx]\n",
    "    \n",
    "    cloned_clf.fit(boot_train_X, boot_train_y) \n",
    "    train_acc = clf.score(X_train, y_train)\n",
    "    valid_acc = cloned_clf.score(boot_valid_X, boot_valid_y)\n",
    "    \n",
    "    gamma = (no_information_rate(\n",
    "        y, cloned_clf.predict(X),\n",
    "        accuracy_score))\n",
    "    R = (valid_acc - train_acc) / (\n",
    "        gamma - train_acc)\n",
    "\n",
    "    weight = 0.632 / (1 - 0.368*R)\n",
    "    \n",
    "    \n",
    "    acc = (weight*train_acc + (1. - weight)*valid_acc)\n",
    "    \n",
    "    bootstrap_train_accuracies.append(acc)\n",
    "\n",
    "bootstrap_train_mean = np.mean(bootstrap_train_accuracies)\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dc273eb5-62bc-40f7-853d-c77346ed4840",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\nci_lower = np.percentile(bootstrap_train_accuracies, 2.5)\\nci_upper = np.percentile(bootstrap_train_accuracies, 97.5)\\n\\n\\nresults['Bootstrap, .632'] =            {'Test accuracy': bootstrap_train_mean,\\n            'Lower 95% CI': ci_lower,\\n            'Upper 95% CI': ci_upper}\\n\\n\""
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "ci_lower = np.percentile(bootstrap_train_accuracies, 2.5)\n",
    "ci_upper = np.percentile(bootstrap_train_accuracies, 97.5)\n",
    "\n",
    "\n",
    "results['Bootstrap, .632'] = \\\n",
    "           {'Test accuracy': bootstrap_train_mean,\n",
    "            'Lower 95% CI': ci_lower,\n",
    "            'Upper 95% CI': ci_upper}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19cc7408-461a-4e39-90ee-5cd7bf48c8af",
   "metadata": {},
   "source": [
    "### 3) Bootstrapping the Test Set predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f1ef995a-65c3-402d-8664-db29cb67f576",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X_train, y_train)\n",
    "\n",
    "predictions_test = clf.predict(X_test)\n",
    "acc_test = np.mean(predictions_test == y_test)\n",
    "\n",
    "rng = np.random.RandomState(seed=12345)\n",
    "idx = np.arange(y_test.shape[0])\n",
    "\n",
    "test_accuracies = []\n",
    "\n",
    "for i in range(200):\n",
    "\n",
    "    pred_idx = rng.choice(idx, size=idx.shape[0], replace=True)\n",
    "    acc_test_boot = np.mean(predictions_test[pred_idx] == y_test[pred_idx])\n",
    "    test_accuracies.append(acc_test_boot)\n",
    "\n",
    "bootstrap_train_mean = np.mean(test_accuracies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6f6fdfb8-38d0-45a6-a554-ee25543e2d4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ci_lower_3 = np.percentile(test_accuracies, 2.5)\n",
    "ci_upper_3 = np.percentile(test_accuracies, 97.5)\n",
    "\n",
    "results[\"Method 3: Bootstrap test set\"] = {\n",
    "    \"Test accuracy\": bootstrap_train_mean,\n",
    "    \"Lower 95% CI\": ci_lower_3,\n",
    "    \"Upper 95% CI\": ci_upper_3,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb4ec325-61b9-4862-acf4-0bfa523d228d",
   "metadata": {},
   "source": [
    "## Comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3be3ae62-a57d-4183-bcf4-b80593e1f7e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAADQCAYAAAD4dzNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzeklEQVR4nO3deZxUxbn/8c9XQBkWxYXwMhjFYBxjgAyRuASVwWvEFUFNcEkUl6hJNDHGBWMu4WISMHoTl1y3+EOMcSdIUDS4wKBBBWWbcZu44YImSHBUYGQZnt8fVa2Hpnu6Z2PmwPN+veZld51zqp6qbuc5VecwR2aGc84559Jlq9YOwDnnnHMN5wncOeecSyFP4M4551wKeQJ3zjnnUsgTuHPOOZdCnsCdc865FPIE7pzbLEkaKOlVSSskDZP0iKTT8uzbS5JJar+p48wnxv3l1o6jEEnlkt5t7Ti2RJ7AnXMASPqqpBmSPpL0mqThiW2ZBLci8fPfie0nS3pf0puSyhPlvSU9Landpu0NAGOBP5pZFzObYmZHmNntrRBHo8S43yhm3/jZ7NHSMbm2pc2cbTrnWk+cef4NuAn4NjAIeFBSfzP7Z2LXbma2Lsex44FvAPsAfwT6xM3XAReaWV0LdyGX3YAXW6HdVJHUPvszdengM3DnHMBewBeBP5hZnZnNAGYD3y/i2B2BJWb2PvA48GUASSfE8mcLVSDpB5JelvSJpJckfSOWf1VShaQaSS9KGpo4ZqKk/5M0LR43R1LvuO31GMeDcbVgm1jPWXF7O0lXS1om6Q3gqKx4tpP0/+KqwhJJv86sIkgaKekf8fgP46rDEYljd5B0m6T34vYpiW1HS1oY+/O0pH71jMlns+oCfX0yHrIo9nVEobYkLZZ0qaRKYKWkX0qalNX+tZKui69PT3w+b0g6p564L41j9omkakn/lW9f10Rm5j/+4z9b+A/QF1gBKFH2GPBAfN0LMGAJ8C5wG7BT3LYV8E9gF+AY4DmgC7AQ2LGItr8T6/0mIGAPwuy5A/Aa8Atga+AQ4BOgNB43EVgO7EtYTbwTuCdR72Lg0MT7CuCs+Ppc4BXgS8AOwMzYv/Zx+xTgZqAz8AVgLnBO3DYSWAv8AGgH/BB4LzN2wDTgXmD72IdBsfwbwFJgv3jcaTHGbfKMiwF7FNnXz/Ytpq34emHsf0kc71XAtnF7O+B9YP/4/iigd/x8BsV9vxG3lQPvxtelwDvAFxPfm96t/f3eXH98Bu6cg5DMlgIXS+og6TDCL+pOcfsyQoLdjbBM3pWQRDCz9YQkNgm4iJDYxgLXA30lzZQ0XVIfcjsL+J2ZPWfBa2b2FrA/4URgvJmtsbAq8BBwUuLYyWY218IS8J1AWZH9/S5wjZm9Y2bLgXGZDZJ6AEcAF5jZSjNbCvwBODFx/Ftm9icLlwZuB3YGekjaOR57rpl9aGZrzWxWPOYHwM1mNsfCKsftwOrYz2I0pK/FtHVd7H9tHO/5wLC47RBglcXVEzObZmavx89nFvAocFCOduuAbYC9JXUws8Vm9nqR/XMN5AncOYeZrSX88j4K+Bfwc+A+wmwbM1thZs+b2Toz+zdwHnCYpG3j9ifMbH8zGwSsBwYQZo13EGasVwC35mn+S0CuX/JfBN6JJwgZbwE9E+//lXi9ipDwi/FFwkwxWW9GZvb/flx+riHMxr+Qq10zWxVfdiH0ZbmZfZijzd2An2fqjPV+KcZSjIb0tZi23sk65i4+Pzk6Ob4HQNIRkp6VtDzWdSSwU3ajZvYacAEwBlgq6R5JxfbPNZAncOccAGZWaWaDzGxHMxtCuIY8N9/u8b9KFkoS4Sa2nxB+wbeLs7vngHzXe98hLM9mew/4kqTk76ldCcvtTfU+IaEl603Gs5pwiaBb/NnWzL5WRL3vADtI6pZn228SdXYzs05mdndjO1EgjkJtZT+K8n6gXNIuwHBiApe0DfBX4Gqgh5l1Ax4m67P/rFKzu8zsQMJJhAFXNmO/XIIncOccAJL6SeooqZOkiwjLwhPjtv0klUraStKOhLvLK8zso6xqzgIWmNlC4D9AiaS9gcFAvn8SdStwkaR9FOwhaTdgDrASuCQu65cTrrHf0wzdvQ/4iaRdJG0PjMpssHAz3qPA/0raNva5t6RBhSqNxz4C3CBp+xj3wXHzn4Bz41hKUmdJR0nq2gz9+Tfx5sHGtmVmHxDuE7gNeNPMXo6btiYsi38ArIs37B2Wq474HTkkJv1PgVrCsrprAZ7AnXMZ3yfMTJcC/wV828xWx21fBv5OuInsBcIMNXktGkk7AT8F/hsgXqs9D5hB+Odp5+dq1MzuB35DmPF9QriBbAczWwMMJVxTXgbcAJxqZq80Q1//BEwHFhGu/U7O2n4qIXG9BHxIuL6/c5F1f59wk1vmvoILAMzsecK16T/GOl8jXF5oDmOA2+Ny+Xeb0NZdwKEkls/N7BPCisp9sa6Tgal5jt+G8E8KlxGW/L9AuAnRtYDMXZPOOeecSxGfgTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUsgfZuK2WN26dbM99kjvA5xWrlxJ586dWzuMRvP4W1ea409z7ADz5s1bZmbdm1qPJ3C3xerRowfPP/98a4fRaBUVFZSXl7d2GI3m8beuNMef5tgBJL1VeK/CfAndOeecSyFP4M4551wKeQJ3zjnnUsgTuHPOOZdCnsCdc865FPIE7pxzxSgvDz/OtRGewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRTyBO6cc86lkCdw55xzLoU8gTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUsgTuHPOOZdCnsCdc865FPIE7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSKFUJXJJJuiPxvr2kDyQ9VOC4MklHJt6PkXRRE+LIebykcyVVSVoo6R+S9i6irgpJ1fGYlyWd3YS4flHEPhdI6tTYNgrUvcE4N/DYbpJ+1ExxjJT0xeaoy7lcpixYws8rVrH7qGkMHD+DKQuWtHZIbguUqgQOrAT6SCqJ778NFPN/ThnQqMTSQHeZWV8zKwN+B/y+yONOiccMBK6UtHUj2y+YwIELgJwJXFK7RrabUUbjx7kb0CwJHBgJeAJ3LWLKgiVcNrmK/3xqGLCkppbLJld5EnebXPvWDqARHgGOAiYBJwF3AwcBSOoMXA/0JfRtTNx/LFAi6UBgXKxnb0kVwK7ANWZ2XazjQuCMuM+tZnZNLL8cOBV4B/gAmJcdmJl9nHjbGbAG9q0L4SSlLrZ5EiEpC5hmZpfmK5c0PvZxIfAicDZwH7AL0A64AuhBSGwzJS0zs8GSVhBONIYAP5d0CHAMUAI8DZxjZhbHaiGwL7AtcIaZzc0EHk86ssf5IbI+DzP7m6SvAbcBWxNOIo+P8fWO8T9mZhcn6u6c3Rczu1fSPjH2LsAyQuIeCAwA7pRUCxxgZrUN/BzcZmbEzc80uY7R74X/vS+ZVMmauvUbbKtdW8clkyq5e+7bTW7n3nMOaHIdbsuQxgR+DzA6Lpv3AyYQEzhwOTDDzM6Q1A2YCzwOjAYGmNl5EJbAgb2AwUBXoFrSjbG+04H9CMlxjqRZhCRzItCfMGbzyZHAY90/Bi4kJKdDEuUL4yw7lzslrQa+AlxgZnVxCfhKYB/gQ+BRScNinzYqN7NRks7LtCHpeOA9Mzsqvt/OzD6KJyiDzWxZbLsz8IKZjY77vWRmY+PrO4CjgQcz+5rZtyQdHMe9T6YDZrZGUvY4/zb785D0OHAucK2Z3RkTfztgFNAnzxgdnt0XSR0IJwfHmtkHkkYAv4ltnQdcZGbP5/h8ziac3NC9e3cqKiryfCRt34oVKzz+ItXUNP0cbt26dQAbJe+MNXXrqampaXI7m2pM0vz9SXPszSl1CdzMKiX1Isy+H87afBgwNHF9uiNhhp3LNDNbDayWtJQwOz0QeMDMVgJImkw4Odgqlq+K5VPrie//gP+TdDLwS+C0WF5WT7dOMbPnJXUHnpb0d8JydIWZfRDbvBM4mDCrz1U+JavOKuBqSVcCD5nZU3nargP+mng/WNIlhGX2HQiz+UwCvzv25UlJ20rqZmY19fQr3+fxDHC5pF2AyWb2qqR6qtm4L5L6EE4gHovHtgPer6+SGPstwC0ApaWlVl5eXuiQNquiogKPvzjN0swjVwLQs1sJS3KcEPTsVsL0Sw/ZqLytSvP3J82xN6e0XQPPmApcTUwoCQKON7Oy+LOrmb2cp47Vidd1hJOZ+rJIQ5fD7wGGNeSAmJTn8/kKQC71ZrpEXf8kzNKrgHFxdpzLp2aWWbLvCNwAnGBmfYE/EZLuZ9VmN1MgjJyfh5ndBQwFaoHpcdm+oX0R8GKi7r5mdliBeJxrsouHlFLSYcPbRUo6tOPiIaWtFJHbUqU1gU8AxppZVVb5dOB8xSmZpP6x/BPCUnkhTwLDJHWK112HA0/F8uGSSiR1JVwj3oikryTeHgW8WmyH4vGdCMv0rwNzgEGSdoo3l50EzKqnHGBtXFomLsGvMrO/EE52vhH3qW8sMsl6maQuwAlZ20fEug8EPjKzj7K2Z9ed8/OQ9GXgjXjfwVTCpYu8ceXpSzXQXdIBcZ8O8dp6oT461yTD+vdk3HF92bGjEGHmPe64vgzr37O1Q3NbmNQtoQOY2bvAtTk2XQFcA1TGpLGYcA13JjAq3iA1LsdxmXrnS5pIuM4M4Sa2BQCS7iXcxPUWIanncp6kQ4G1hOvTp2U2FHENvBbYBphoZvPiMZfF2AU8bGZ/q6+csDRcKWk+8GfgKknrYzw/TOzziKT3zWxwVv9rJP2JMNNdDDyXFeeHkp4m3sSWox/Z45zv8xgBfE/SWuBfhJOx5ZJmS3oBeCR5ExvhJrgN+hKvuZ8AXCdpO8J3+RrCkv9E4Ca/ic21lGH9e9Lto1d9Gde1Kpk1dGXYbYniXeg5bwxLq9LSUquurm7tMBot7dcBUxd/JtZ481Tq4s+S5vjTHDuApHlmNqCp9aR1Cd0555zboqVyCd1temZW3toxOOec+5zPwJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRTyBO6cc86lkCdw55xzLoU8gTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUiiVCVySSboj8b69pA8kPVTguDJJRybej5F0URPiyHm8pAslvSSpUtITknbLc/zfJS2S9KKkmyS1K9BeL0m1khbG456WVNrI2DcYi3raO7kx9TeVpH6SnoljUyWpYyzPOWbFjrlzzm0uUpnAgZVAH0kl8f23gSVFHFcG1Ju0mskCYICZ9QMmAb/Ls993zezrQB+gO/CdIup+3czK4nG3A79oZIxlFB6LXkDOBC6pfSPbLSjW/RfgXDP7GlAOrI2b841ZsWPuHABTFixh4PgZ7D5qGgPHz2DKgmJ+hTjXdqQ1gQM8AhwVX58E3J3ZIKmzpAmSnpO0QNKxkrYGxgIj4gx2RNx9b0kVkt6Q9JNEHRdKeiH+XJAov1xStaTHgZyzXzObaWar4ttngV3y7PdxfNke2BqwBo7BtsCHMa6Okm6Ls9UFkgbnK881FpIGxdcL435dgfHAQbHsZ5JGSrpf0oPAo5K6xNnu/Fj/sbHNXpJekXR7nBFPktSpAf06DKg0s0VxnP5jZnX1jVmxY+4chOR92eQqltTUYsCSmloum1zlSdylSovNojaBe4DRcdm8HzABOChuuxyYYWZnSOoGzAUeB0YTZmnnQVgCB/YCBgNdgWpJN8b6Tgf2AwTMkTSLcMJzItCfMHbzgXkF4jyTcLKRk6TpwL5xn0mxbGiMc3SOQ3pLWhjj7RRjBPgxgJn1lbQXIcHumasc2DPHWDwI/NjMZkvqAnwKjAIuMrOj4z4jgQOAfma2PM6Uh5vZx5J2Ap6VNDXGUwqcGeubAPwIuLrAWGXsCVgcm+7APWb22Yw615hlqXfM3eZlxM3PNPiYBW/XsKZu/QZltWvruGRSJXfPfTvnMaPfC+eOY2N7NTW1lJc3uGnnmk1qE7iZVUrqRZh9P5y1+TBgaOL6dEdg1zxVTTOz1cBqSUuBHsCBwANmthJA0mTCycFWsXxVLJ+ap07i9u8BA4BB9fRjSLy+eydwCPCYmU0F8tX9upmVxfpHALcAh8eYr491viLpLUIizFeebTbwe0l3ApPN7F1Judp/zMyWZ7oI/FbSwcB6oCdh/ADeMbPZ8fVfgJ9QfAJvH+P+JrAKeELSPDN7IvZjozHLHFhozCWdDZwN0L17dyoqKooMqe1ZsWKFx09IpA2VnbyT5TU1NTm3rVu3LrYXttfV1fn4t5I0x96cUpvAo6mEpFAO7JgoF3C8mVUnd5a0HxtbnXhdRxiTnJkrKmqZW9KhhJWAQfEEIX+FZp/Gk4FjSSSjIkwFbss0mS+UYioys/GSphGuiz8b489lZeL1KYQZ8j5mtlbSYsLJEmw8Tg25PPAuMMvMlgFIehj4BvBEIt6NxqyYMTezWwgnPZSWllp5iqdQFRUVePw0ahY8cPwMluRI/D27lTD90kNyH/TIlQBMv/QIwMe/NaU59uaU5mvgEJbNx5pZVVb5dOB8xSmkpP6x/BPC0nMhTwLDJHWS1BkYDjwVy4dLKonXiI/JdXBs72ZgqJktzbNPF0k7x9ftCYnzlSJiSzoQeD0R8ymxvj0JKw7V9ZRvMBaSeptZlZldCTxPuLRQaLy2A5bG5D0YSN75vaukA+Lrk4B/xHbGSRpeoF/TgX5x/NsTZtMv1TdmxYy5cxkXDymlpMOG/+ijpEM7Lh7SqH/U4VyrSPUM3MzeBa7NsekK4BqgMibxxcDRwExgVLyGPK6eeudLmki4dg5wq5ktAJB0L7AQeIuQ1HO5CugC3B/PId42s6Hx+IVxCbwzMFXSNkA7YAZwU9ynmGvgAtYAZ8XyG4CbJFUB64CRZrZaUr7y7LE4MCbhOuAlwjXk9cA6SYuAicQb5hLuBB6U9Hwck+QJyMvAaZJuBl4FbozlfclxeSDZZzP7UNLvgecIM/eHzWyapB75xqy+MXcu27D+PQG4ano179XU8sVuJVw8pPSzcufSQGYNvfHZufrFexMeMrM+ObZNN7Mhmz6qjZWWllp1dXXhHduotC8jpi7+TKzx2mvq4s+S5vjTHDtAvKdnQFPrSfsSukuZtpK8nXMu7VK9hO7aJjNbTPhDK84551qIz8Cdc865FPIE7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQs2SwCWZpDsS79tL+kDSQwWOK5N0ZOL9GEkXNSGOnMdLulDSS5IqJT0habcc+3SSNE3SK5JelDS+iPZ6SaqVtFDSIklPSyptZOwbjEU97Z3cmPrTKsd3ZKikUfF1k74vzjmXZs01A18J9JFUEt9/G1hSxHFlQL1Jq5ksAAaYWT9gEvC7PPtdbWZ7Af2BgZKOKKLu182szMy+DtwO/KKRMZZReCx6ATkTuKT2jWy32Ulq14zVlZEYFzObamYFT66cy5iyYAkDx89g91HTGDh+BlMWFPOrybm2rzmX0B8BjoqvTwLuzmyQ1FnSBEnPSVog6VhJWwNjgRFxBjsi7r63pApJb0j6SaKOCyW9EH8uSJRfLqla0uNAztmvmc00s1Xx7bPALjn2WWVmM+PrNcD8XPsVsC3wYYyro6TbJFXFPg/OV55rLCQNiq8Xxv26AuOBg2LZzySNlHS/pAeBRyV1iSsM82P9x8Y2e8WVhdvjKsQkSZ3q64ikiZJukvSUpH9KOjqWt5N0VfwsKyWdE8vLJc2UdBdQFfe7OsZRKen8uN8+kmZJmidpuqSdY3mFpCslzY3tHZRnXEZK+mOOeHtL+nus9ylJezXws3OboSkLlnDZ5CqW1NRiwJKaWi6bXOVJ3G0WmnPWdg8wWmHZvB8wATgobrscmGFmZ0jqBswFHgdGE2bG50FYEgX2AgYDXYFqSTfG+k4H9gMEzJE0i3ACciJhxtyekHTnFYjzTMLJRl4xxmOAa+P7oTHO0Tl27y1pYYy3U4wR4McAZtY3JpNHJe2ZqxzYM8dYPAj82MxmS+oCfAqMAi4ys0wyHQkcAPQzs+VxFj7czD6WtBPwrKSpMZ5S4MxY3wTgR8DVBcaqFzAI6A3MlLQHcCrwkZl9U9I2wGxJj8b99wX6mNmbkn4I7A70N7N1knaQ1AG4HjjWzD6IJ22/Ac6Ix7c3s33jkvmvzOxQSdnjMjJPrLcA55rZq5L2A24ADinQP9dCRtz8TMF9ampqubG68H5NseDtGtbUrd+grHZtHZdMquTuuW83qK7R730MwNjYt4bGf+85BzSoPecKabYEbmaVknoRZt8PZ20+DBiqz69XdgR2zVPVNDNbDayWtBToARwIPGBmKwEkTSacHGwVy1fF8ql56iRu/x4wgJCU8u3TnrB6cJ2ZvRH7NhXIV/frZlYWjx1BSCSHx5ivj8e/IuktQqLOV55tNvB7SXcCk83sXUm52n/MzJZnwgd+K+lgYD3QkzB+AO+Y2ez4+i/ATyicwO8zs/XAq5LeIJxcHQb0k3RC3Gc74CvAGmCumb0Zyw8FbjKzdbGvyyX1AfoAj8W+tAPeT7Q3Of53HuHkoSjxBOdbwP2JMdomz75nA2cDdO/enYqKimKbaXNWrFjRZuOvqaktuE9dXR01NTUtGkd28k6WN7TtdevWAXx2XEPjb2ufVVv+/hSS5tibU3NfN51KSArlwI6JcgHHm1l1cuc4U8q2OvG6LsaYM3NFVkxgkg4lrAQMiicI+dwCvGpm1xRTb5apwG2ZJvOFUkxFZjZe0jTC9d9nY/y5rEy8PgXoDuxjZmslLSacLMHG41TMuOU6RsD5ZjY9uUFSeVYsynG8gBfNLN9UJPO5ZD73Ym0F1GROpOpjZrcQPmNKS0utvLy8Ac20LRUVFbTV+IsJa1PEP3D8DJbkOJno2a2E6Zc2cIHmkSsBmH5puDWmLY9/MdIcf5pjb07N/c/IJgBjzawqq3w6cL7i9EhS/1j+CWHpuZAngWEKd4p3BoYDT8Xy4ZJK4jXiY3IdHNu7GRhqZkvzNSLp14QZ5QVFxJTLgcDriZhPifXuSVhxqK6nfIOxkNTbzKrM7ErgecLst9B4bQcsjcl7MJC8235XSZnEeRLwj9jOOEnD89T3HUlbSeoNfDnGOR34YVwOR9Ke8TPJ9ihwblzRQNIO8fjumTgkdZD0tXr6QxF9xsw+Bt6U9J1YryR9vUC9bgtw8ZBSSjpseE9lSYd2XDykUf9YxLk2pVkTuJm9a2bX5th0BdABqJT0QnwPMJNw01ryJrZc9c4HJhKunc8BbjWzBbH8XmAh8FdCUs/lKqALYYl1YXKpPV6/RtIuhBn63sD8uN9ZcdtQSWPz1N077rsI+C1wViy/AWgnqSrGODLO/POVZ4/FBQo37C0CagnX7SuBdQr/ZO1nOWK5Exgg6XnCScIriW0vA6dJqgR2AG6M5X2Bf+XpWzUwK7Z9rpl9CtwKvBTH6AXCiVGu2fKtwNuEz3wRcHK8OfAE4MpYtpCw9F2for4jsb9nxnpfBI4tUK/bAgzr35Nxx/WlZ7cSRJh5jzuuL8P692zt0JxrMpkVtQLtUizem/CQmfXJsW26mQ3JUT4xHjOp5SNsHaWlpVZdXV14xzYq7cuIqYs/E2u89pq6+LOkOf40xw4gaZ6ZDWhqPf6X2LZwuZK3c865tq/N/PEP13LMbDHh7u+GHDOyRYJxzjnXLHwG7pxzzqWQJ3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRQqmMAlmaQ7Eu/bS/pA0kMFjiuTdGTi/RhJFzU20HzHS7pQ0kuSKiU9IWm3PMf/RtI7klYU2V4vSbWSFkpaJOlpSaWNjH2DsainvZMbU38j4jlY0nxJ6ySdsCnaLERSeaHvVI5j9pT0sKTXJL0s6T5JPRpTl9t8TVmwhIHjZ7D7qGkMHD+DKQuWtHZIzjWLYmbgK4E+kkri+28DxfwfUAbUm7SayQJggJn1AyYBv8uz34PAvg2s+3UzKzOzrwO3A79oZIxlFB6LXkDOBC6pfSPbzedtYCRwVzPXu8lI6ghMA240sz3M7KvAjUD31o3MtSVTFizhsslVLKmpxYAlNbVcNrnKk7jbLBSbGB4BjiIkyJOAu4GDACR1Bq4H+sb6xsT9xwIlkg4ExsV69pZUAewKXGNm18U6LgTOiPvcambXxPLLgVOBd4APgHnZgZnZzMTbZ4Hv5eqAmT0b6yyyyxvZFvgw1tGRkCwGAOuAC81sZq5yYDYbj8W/gGszoQEHA+OBr0paSDhZ+JAw5h2BzpKGAn8Dtgc6AL80s79J6gX8HZgD9Af+CZxqZqvydcTMFsd+rK+vw5K+A/wKqAM+MrODY3t3AJ3jbueZ2dOSyoH/Af5NOGGZDFQBPwVKgGFm9rqkicCnwNeAHnHsNpgt5/pOmdnfssI7GXjGzB5M9GtmPL68vn65TWPEzc/Uu72mppYbq+vfp6kWvF3DmroNv+a1a+u4ZFIld899u0F1jX7vYwDGxn41R/z3nnNAk453W7ZiE/g9wOi4LNkPmEBM4MDlwAwzO0NSN2Au8DgwmjAzPg/CEjiwFzAY6ApUS7ox1nc6sB8gYI6kWYTVgRMJSak9MJ8cCTzLmYSTh6LFxDjAzEbn2Nw7JtSuQKcYI8CPAcysr6S9gEcl7ZmrHNgzx1g8CPzYzGZL6kJIaKOAi8zs6LjPSOAAoJ+ZLY+z8OFm9rGknYBnJU2N8ZQCZ8b6JgA/Aq5uyDjkMRoYYmZL4mcLsBT4tpl9KukrhJO5AXHb14GvAsuBNwgnY/tK+ilwPnBB3K8XMAjoDcyUtEdWuxt9pyQ9bmYrE/v0ofD3YSOSzgbOBujevTsVFRUNraLNWLFiRZuOv6amtt7tdXV11NTUtGgM2ck7Wd7QttetWwfw2XHNEX9rfn5t/ftTnzTH3pyKSuBmVhlnXicBD2dtPgwYmrg+3ZEww85lmpmtBlZLWkqYgR0IPJD55SxpMuHkYKtYviqWT81TJ3H79wiJZFAxfUr0bSqQr+7Xzaws1j8CuAU4PMZ8fTz+FUlvERJ1vvJss4HfS7oTmGxm7+ZZGXjMzJZnugj8VtLBwHqgJ2H8AN4xs9nx9V+An9A8CXw2MFHSfYQZNYTZ/x8llRFm5sn+PWdm7wNIep1wAgNhJj44sd99ZrYeeFXSG4QTu6R836mXm9ohM7uF8DlSWlpq5eXlTa2y1VRUVNCW4y8U2qaIf+D4GSzJcSLRs1sJ0y89pGGVPXIlANMvPQJo++NfSJrjT3Pszakhd6FPJSSFu7PKBRwfrxWXmdmuZpbvF+3qxOs6wglEfWvaVkxgkg4lzNqGxhOEljCVsNQN+WMuan3ezMYDZxGWlp+Ns/VckjPOUwjXd/eJJxX/JiQ22Hicihq3bPFGv4Vx1QEzOxf4JfAlYKGkHYGfxba/Tjhh2jpRRXLs1yfer2fDk8VC8RbznXoR2Kch/XNbnouHlFLSod0GZSUd2nHxkEbdj+pcm9KQBD4BGGtmVVnl04HzFaeQkvrH8k8IS8+FPAkMk9QpXvscDjwVy4dLKpHUFTgm18GxvZsJyXtpA/rTUAcCrydiPiW2vydhdlhdT/kGYyGpt5lVmdmVwPOEGWih8doOWGpmayUNBpJ32+8qKXMx7STgH7GdcZKGF9tBM7s8kzQTcc6JlxeWERL5dsD7cQb9faBd3grz+46krST1Br5MGKOkfN+ppLuAb0k6KlMg6XBJfRsRj9tMDevfk3HH9aVntxJEmHmPO64vw/r3bO3QnGuyou9uNrN3+fzGq6QrgGuAyvgLdzFwNDATGBVnc+NyHJepd368sWluLLrVzBYASLoXWAi8RUjquVwFdAHuj7/v3zazofH4hYlk9DvCjU+dJL0b2xlT5DVwAWsIs2aAG4CbJFURblYbaWarJeUrzx6LA2MSrgNeIly3Xw+sk7QImEi8YS7hTuBBSc/HMXklse1l4DRJNwOvEm6kg3AT2EaXByR9E3iAcEPcMZL+x8y+lmts43VuAU8Ai2Lf/xpvcJvJhqsExaoGZhEuAZwbr6cnt+f7Tn3GzGolHQ1cI+kaYC1QSbhpbsdGxOQ2U8P69/SE7TZLMmvUaqtrI+K9CQ+ZWZ8c26ab2ZBNH1V+8WTtITOb1NqxlJaWWnV19uQ/PdJ+HTB18WdijTdPpS7+LGmOP82xA0iaZ2YDCu9ZP/9LbJuxtpa8nXPONZ/m/gMhbhOL/6Z7o9l3W2VmI1s7Buec2xz4DNw555xLIU/gzjnnXAp5AnfOOedSyBO4c845l0KewJ1zzrkU8gTunHPOpZAncOeccy6FPIE755xzKeQJ3DnnnEshT+DOOedcCnkCd84551LIE7hzzjmXQp7AnXPOuRTyBO6cc86lkCdw55xzLoU8gTvnnHMp1OwJXJJJuiPxvr2kDyQ9VOC4MklHJt6PkXRRE+LIebykgyXNl7RO0glF1lUh6fnE+wGSKhobW2M0dTxamqShkkY1U10XSOqUeP+wpG7NUbdzzWHKgiX8vGIVu4+axsDxM5iyYElrh+S2QC0xA18J9JFUEt9/Gyjm210GHFlop2bwNjASuKuBx31B0hGNaVBS+8Yc15oktWvI/mY21czGN1PzFwCfJXAzO9LMapqpbueaZMqCJVw2uYr/fGoYsKSmlssmV3kSd5tcSyWWR4CjgEnAScDdwEEAkjoD1wN9Y/tj4v5jgRJJBwLjYj17x5nursA1ZnZdrONC4Iy4z61mdk0svxw4FXgH+ACYlx2YmS2O+65vYJ+uAn4ZY/2MpI7AjcAAYB1woZnNlDQyjkFHoLOkPwPDgHZAH+B/ga2B7wOrgSPNbLmkHwBnx22vAd83s1X5gpJ0TIxra+A/wClm9m9JY4DeQE/gS8DvzOxPksoJY/0foBR4EviRma2XtAL4PTAE+Lmkfcka5zj2fczsDEl9CZ/tvsB3gQFmdp6kiUAtsBewG3A6cBpwADDHzEbG2G8EvgmUAJPM7FeSfgJ8EZgpaZmZDZa0ONa9LNdnL6lX/Fz+AXyLcMJ4rJnV5hs3t2UZcfMzTa5j9HsfA3DJpErW1G3466N2bR2XTKrk7rlvN7mde885oMl1uC1DSyXwe4DRcdm8HzCBmMCBy4EZMQF0A+YCjwOjiQkAwpIxIQEMBroC1fEXfj9CQtgPEDBH0izCasKJQP/Yr/nkSOD1kbTQzMrybH4GGC5pMPBJovzHAGbWV9JewKOS9ozbDgD6xcQ8kpC4+xOS+mvApWbWX9IfCCce1wCTzexPMZ5fA2cSTnjy+Qewv5mZpLOAS4Cfx239gP2BzsACSdNi+b7A3sBbwN+B4wgnW52BF8xstKR9yD3O1wAVkoYTPstzzGyVpOy4tgcOAYYCDwIDgbOA5ySVmdlC4PI4Nu2AJyT1M7PrYpIebGbLkhXWE9OHwFeAk8zsB5LuA44H/pIdlKSzCSdIdO/enYqKinqGtm1bsWKFx1+kmpqmn8tdeOZvAFjzYe5z/zV166mpqWlyO5tqTNL8/Ulz7M2pRRK4mVXGWdFJwMNZmw8Dhiau53YkzLBzmWZmq4HVkpYCPYADgQfMbCWApMmEk4OtYvmqWD61EXGXFdjl14TZ7qWJsgOJCdbMXpH0FpBJ4I+Z2fLEvjPN7BPgE0kfERIbQBUh2UK4/PBroBvQBZheIKZdgHsl7UyYhb+Z2Pa3OAutlTSTkLhrgLlm9gaApLtjHyYBdcBfE/3aaJzNbEE8GakEbjaz2XniejCeVFQB/zazqljPi0AvYCHw3ZhQ2wM7E04qKuvpa77PfirwZjwpgHDi1itXBWZ2C3ALQGlpqZWXl9fTXNtWUVGBx1+c5mxm4PgZLMlxQtCzWwnTLz2k+RpqYWn+/qQ59ubUknehTwWuJiyxJgk43szK4s+uZvZynjpWJ17XEX7RbzTVS7BGR1sEM5tBOOHYP1FcXzwrs94n+7M+8X49n59MTQTOM7O+wP/E9upzPfDHuP85Wftnj4cVKP/UzOri6/r69RVgBWGpO59k37L73V7S7sBFwH+ZWT9gGoX7Wl9Mub4rzjW7i4eUUtJhw1tESjq04+Ihpa0UkdtStWQCnwCMzcy8EqYD5yuuuUrqH8s/ISyVF/IkMExSp3g9fTjwVCwfLqlEUlfgmOboRA6/ISxTJ+M5BSAune8KVDeh/q7A+5I6ZOotYDs+v0nwtKxtx0rqKGlHoBx4LpbvK2l3SVsBIwjL8NlyjrOk7YBrgYOBHYu9kz+HbQknOB9J6gEkbxDM913I99k7t8kM69+Tccf1ZceOQoSZ97jj+jKsf8/WDs1tYVpslmJm7xJ+0We7gnAdtTIm8cXA0cBMYJSkhXx+E1uueufHm6TmxqJbzWwBgKR7CUuzb5HnF7ukbwIPEK7RHiPpf8zsa3FbfdfAM+0/LOmDRNENwE1xqXgdMNLMVue4Jlys/wbmxD5UUfikZgxwv6QlwLPA7oltcwkz212BK8zsvXiS8QwwnnAj4ZOE8dhAvnGWNAG4wcz+KelMws1mTza0k2a2SNIC4EXgDSC5FH8L8Iik981scBEx9Wpo+841xbD+Pen20au+jOtalcxadNXZtZJ4E+AKM7s6q7wcuMjMjm6FsNqU0tJSq65uymJJ60r7dUCPv3WlOf40xw4gaZ6ZDWhqPf6X2JxzzrkU8hm422JJ+oSm3a/Q2nYClhXcq+3y+FtXmuNPc+wApWZWzD1f9fI7dd2WrLo5lrFai6TnPf7W4/G3njTHDiH+5qjHl9Cdc865FPIE7pxzzqWQJ3C3JbultQNoIo+/dXn8rSfNsUMzxe83sTnnnHMp5DNw55xzLoU8gTvnnHMp5AncbXYkHS6pWtJrkkbl2H6xpIXx5wVJdZJ2KObYTaGJ8S+WVBW3Ncs/VWmoIuLfTtKDkhZJelHS6cUeuyk0Mf40jP/2kh6QVClprqQ+xR67KTQx/lYdf0kTJC2V9EKe7ZJ0XexbpaRvJLY1fOzNzH/8Z7P5AdoBrwNfJjxedRGwdz37H0N4Pn2Dj21r8cf3i4Gd2vL4A78AroyvuwPL476pGP988ado/K8CfhVf7wU80ZjvXluLv42M/8HAN4AX8mw/EniE8GTF/YE5TRl7n4G7zc2+wGtm9oaZrQHuAY6tZ/+T+PyRtw09tiU0Jf62oJj4Deiq8MSfLoQEuK7IY1taU+JvC4qJf2/gCQAzewXopfBEwLSMf774W52ZPUn4PuRzLPBnC54FuknamUaOvSdwt7npCbyTeP9uLNuIpE7A4cBfG3psC2pK/BCSy6OS5kk6u8WizK+Y+P8IfBV4j/DEvZ+a2foij21pTYkf0jH+i4DjACTtC+wG7FLksS2tKfFD649/Ifn616ix9z+l6jY3uZ7jmu/fSh4DzDazzBlzQ45tKU2JH2CghcfGfgF4TNIrcVawqRQT/xDCY38PAXoT4nyqyGNbWqPjN7OPScf4jweuVXh0cxWwgLCCkJbxzxc/tP74F5Kvf40ae5+Bu83Nu8CXEu93IcyUcjmRDZefG3JsS2lK/JjZe/G/SwnPed+3BWKsTzHxnw5MjsuIrwFvEq5lpmX888WfivE3s4/N7HQzKwNOJVzHf7OYYzeBpsTfFsa/kHz9a9zYt9bFfv/xn5b4IawqvQHszuc3g3wtx37bEa5VdW7osW04/s5A18Trp4HD21r8wI3AmPi6B7CE8HSpVIx/PfGnZfy78flNdz8gXJNNzfe/nvhbffxj273IfxPbUWx4E9vcpoy9L6G7zYqZrZN0HjCdcGfnBDN7UdK5cftNcdfhwKNmtrLQsWmJn5BMHgj3VtEeuMvM/r7poi86/iuAiZKqCL/ILjWzZQApGf+c8Uv6MukY/68Cf5ZUB7wEnFnfsWmJnzbw/Zd0N1AO7CTpXeBXQIdE7A8T7kR/DVhFWM1p9Nj7n1J1zjnnUsivgTvnnHMp5AncOeecSyFP4M4551wKeQJ3zjnnUsgTuHPOOZdCnsCdc61C4SlqmSeq3R//NGxj65oo6YT4+lZJe9ezb7mkbyXenyvp1Ma27Vxr8QTunGsttWZWZmZ9gDXAucmNkto1plIzO8vMXqpnl3LgswRuZjeZ2Z8b01ZrkeR/w8N5AnfOtQlPAXvE2fFMSXcBVZLaSbpK0nPx+cnnwGfPVf6jpJckTQO+kKlIUoWkAfH14ZLmKzy7+wlJvQgnCj+Ls/+DJI2RdFHcv0zSs7GtByRtn6jzyvj86X9KOii7A5K6xDbmKzyT+tjEtlNjnYsk3RHLesQ2FsWfb0nqpcSzpCVdJGlMIobfSpoF/FTSMZLmSFog6XHFJ3LFOG6LMVRKOl7SmZL+kKj3B5J+3zwfnWstfhbnnGtVcTZ5BJD5q1n7An3M7E2FJ0p9ZGbflLQNMFvSo0B/oBToS/gLXC8BE7Lq7Q78CTg41rWDmS2XdBOwwsyujvv9V+KwPwPnm9ksSWMJf0nrgritvZntK+nIWH5oVlc+BYab2ceSdgKelTSV8PjLywkP2lgmaYe4/3XALDMbHlcbugDbFxiubmY2KMa9PbC/mZmks4BLgJ8D/x3HrG9ivzVApaRLzGwt4S+AnVOgLdfGeQJ3zrWWEoUnSkGYgf8/wtL2XDN7M5YfBvTLXN8m/A34rwAHA3ebWR3wnqQZOerfH3gyU5dt+NS2jUjajpAgZ8Wi24H7E7tMjv+dR/h71xtVAfxW0sFA5vGoPQhPLZuU+XOxiTgOITyMg9iPjzIz/nrcm3i9C3CvwvOktyY+0INwYnFiZicz+zD2bwZwtKSXgQ5mVlWgLdfGeQJ3zrWWWgtPlPpM/DvWyb/vLsKMeHrWfkdS+HGLKmKfhlgd/1tH7t+dpxCejLWPma2VtBjo2MA41rHhpc2OWduTY3M98HszmyqpHBgTy/O1dyvwC+AV4LYi43FtmF8Dd861ZdOBH0rqACBpT0mdgSeBE+M18p2BwTmOfQYYJGn3eGxm6foToGv2zmb2EfBh4vr294FZ2fvVYztgaUzeg4HdYvkTwHcl7ZgVxxPAD2NZO0nbAv8GviBpx3jJ4OgC7S2Jr09LlD8KnJd5k5nVm9kcwiMrTybrMbQunTyBO+faslsJ17fnx5u7bibMfh8AXgWqCI/33CjRmtkHwNnAZEmL+Hz5+UFgeOYmtqzDTgOuklQJlAFjGxDrncAASc8TZuOvxDheBH4DzIpxZG4e+ykwWOGpZvMIj49cG9ucAzyUqSOPMcD9kp4CliXKfw1sr/DP8xax4cnNfcDszLK6Szd/Gplzzm0hJD0E/MHMnmjtWFzT+QzcOec2c5K6Sfon4b4DT96bCZ+BO+eccynkM3DnnHMuhTyBO+eccynkCdw555xLIU/gzjnnXAp5AnfOOedS6P8DaGhUE3F7T9kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "labels = list(results.keys())\n",
    "\n",
    "means = np.array([results[k][\"Test accuracy\"] for k in labels])\n",
    "lower_error = np.array([results[k][\"Lower 95% CI\"] for k in labels])\n",
    "upper_error = np.array([results[k][\"Upper 95% CI\"] for k in labels])\n",
    "\n",
    "asymmetric_error = [means - lower_error, upper_error - means]\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 3))\n",
    "ax.errorbar(means, np.arange(len(means)), xerr=asymmetric_error, fmt=\"o\")\n",
    "ax.set_xlim([0.7, 1.0])\n",
    "ax.set_yticks(np.arange(len(means)))\n",
    "ax.set_yticklabels(labels)\n",
    "ax.set_xlabel(\"Prediction accuracy\")\n",
    "ax.set_title(\"95% confidence intervals\")\n",
    "\n",
    "ax.vlines(acc_test_true, [0], 5, lw=1.5, color=\"red\", linestyle=\"-\", label=\"True value\")\n",
    "\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"matplotlib-figures/comparison-simulation.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5011b685-be7d-42ac-b2a1-ca6d5a8a2d37",
   "metadata": {},
   "source": [
    "- As it turns out all methods' 95% confidence intervals contain the true test accuracy, which is great. However, I noticed that changing the parameters of the data generating function slightly will have a huge effect on this outcome, so take it with a grain of salt.\n",
    "\n",
    "- Ok, ideally we also want to repeat this simulation many times and see if the confidence intervals really contain the true parameter 95% of the time. I originally wanted to leave this as an exercise to the reader, but then I couldn't resist 😛; see [ci-simulation-repeated.ipynb](ci-simulation-repeated.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
